Relaxation of Terrace-width Distributions: Physical Information from Fokker-Planck 
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Recently some of us have constructed a Fokker-Planck formalism to describe the equilibration of 
the terrace-width distribution of a vicinal surface from an arbitrary initial configuration. However, 
the meaning of the associated relaxation time, related to the strength of the random noise in the 
underlying Langevin equation, was rather unclear. Here we present a set of careful kinetic Monte 
Carlo simulations that demonstrate convincingly that the time constant shows activated behavior 
with a barrier that has a physically plausible dependence on the energies of the governing microscopic 
model. Furthermore, the Fokker-Planck time at least semiquantitatively tracks the actual physical 
time. 

PACS numbers: 05.10.Gg, 68.35.-p, 81.15.Aa, 05.40.a 
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INTRODUCTION 

With equilibrium properties of vicinal surfaces — 
especially the form of the terrace width distribution 
(TWD) — now relatively well understood much atten- 
tion is focusing on non-equilibrium aspects, which have 
long been of interest. In a previous paper some of us 
derived the following Fokker-Planck (FP) equation 
[Eq. |T])] to describe the distribution of spacings between 
steps on a vicinal surface during relaxation to equilib- 
rium. The goal was to describe the relaxational evolu- 
tion of this spacing distribution rather than the evolu- 
tion of the positions of individual steps as in a previous 
investigation [1, 0, H, @| . As in all those papers, we sim- 
plify to a one-dimensional (ID) model, in which a step is 
represented by its position in the x direction (the down- 
stairs direction in "Maryland" notation), averaged over 
the y direction (along the mean direction of the step, the 
"time-like" direction in fermionic formulations) 0]. This 
picture implicitly assumes that one is investigating time 
scales longer than that of fluctuations along the step. 

We started with the Dyson Coulomb gas/Brownian 
motion model; 0, Q made the mean- field- like assump- 
tion, when computing interactions, that all but adja- 
cent steps are separated by the appropriate integer mul- 
tiple of the mean spacing; and set the width of the con- 
fining [parabolic] potential in the model to produce a 
self-consistent solution. Details are provided in the Ap- 
pendix, which expands the earlier derivation and corrects 
some inconsequential errors in intermediate stages Q. 
We found the following: 



the generalized Wigner surmise (GWS), thus 
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where s is the distance w between adjacent steps divided 
by its average value (w) , determined by the slope of the 
vicinal surface. 

The steady-state solution of Eq. (fT]) has the form of 



P e (s) 



(2) 



2[r(ffi]*- 1 



2 



where the constants b g and a g assure unit mean and nor- 
malization, respectively. (The Wigner surmise, Eq. (2) 
pertains to the special cases 0=1,2,4; the generalization 
is to use this expression for arbitrary g > 1.) The 
dimcnsionlcss variable g gauges the strength A of the 
A/w 2 energetic repulsion between steps: (g — l) 2 = 
1 + AA(3 / (k B T) 2 , where (3 is the step stiffness. The di- 
mensionless FP time t can be written as i/r; here the 
relaxation time r is (w) 2 /T, where y/T is the strength of 
the white noise in the Langevin equation (for the step 
position) underlying the FP equation 

To confront data, both experimental and simulational, 
one typically investigates the variance <r 2 (t) or standard 
deviation o~(t) of this distribution. If the initial configura- 
tion of the vicinal surface is "perfect" (i.e. has uniformly- 
spaced straight steps), then a(t) obeysfH 
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where a 2 at = 



cr(oo) is the variance for an infinite system 
at long time. When dealing with numerical data, we take 
the variance to be normalized by the mean spacing, so 
divided by the squared mean terrace width, to mimic the 
formal analysis. The precise value of the proportional- 
ity constant is not of importance to our analysis, since 
we view r as the source of information for an activated 
process, with any prefactor therefore insignificant. As 
discussed in the appendix (esp. Eq. (|A15[) b one might 
expect the prefactor to be unity when the first moment 
has the assumed GWS value of one, but with the ap- 



2 



proximations we make to obtain a compact solution, the 
prcfactor seems better described as two. 

Time in this formulation is not the natural fermionic 
time associated with the direction along the steps (y 
in "Maryland notation"), i.e., that resulting from the 
standard mapping between a 2D classical model and a 
(1+1)D quantum model. Instead it measures the evo- 
lution of the 2D or (1+1)D system toward equilibrium 
and the thermal fluctuations underlying dynamics. Since 
the time constant r enters rather obliquely through the 
noise force of the Langevin equation, a key investiga- 
tional objective in the previous Letter[2j and in this pa- 
per is whether r corresponds to a physically significant 
rate. Monte Carlo simulations allow the examination of 
a well-controlled numerical experiment. In the former 
we used our well-tested Metropolis algorithm to study a 
terrace-step-kink (TSK) model of the surface. We found 
a satisfactory fit to the form of Eq. (J3j> , from which we 
obtained r « 714 MCS (Monte Carlo steps per site) for 
g = 2 (or A = 0, only entropic repulsions) while r w 222 
MCS for g sa 4.47. This result is in qualitative agreement 
with the understanding that T should increase (and, so,r 
should decrease) with increasing g, as discussed in Ref.H. 

In the present paper, we confront more systematically 
and thoroughly the above-noted crucial issue, showing 
that the time constant associated with the FP transcrip- 
tion can be related to the atomistic processes underlying 
the relaxation to equilibrium and that the FP time in 
some sense tracks (though of course does not replicate) 
the literal physical time of the relaxing system. We use 
a standard, simple lattice model that embodies the ba- 
sic atomistic properties of these surfaces. We report far 
more extensive simulations, using kinetic Monte Carlo 



(KMC)[10|, 111] rather than the Metropolis algorithm, for 



a solid-on-solid (SOS) rather than a TSK model, so that 
we have real mass transport. Since atomic energies in 
this generic model are proportional to the number of lat- 
eral nearest neighbors, detailed-balance is satisfied. To 
simplify the analytic expressions and, especially, the sim- 
ulations, we concentrate in this paper on the special case 
g = 2, corresponding to steps with only entropic repul- 
sions ("free fermions"). Wc find that the time constant, 
extracted from the numerical data by fitting to the time 
correlation function in the form predicted by the FP anal- 
ysis, has an activated form that can be related to an 
atomistic rate-limiting process in the simulations. Our 
goal is not to find the best accounting for the dynamics 
of a real stepped surface, nor even of our model surface. 
It is to show that the FP approach offers a relatively 
simple and physically viable approach to accounting for 
the relaxation of artificial initial configurations toward 
equilibrium. 

The second section describes the model and KMC al- 
gorithm that we use. The third presents our numerical 
results. The fourth discusses them, with one subsection 
describing the crucial role played by the creation of kink- 



antikink pairs and another investigating the evolution of 
the shape of the distribution. The fifth makes compar- 
isons with the venerable mean-field treatment of step dis- 
tributions, and the final section sums up our findings. In 
an Appendix we expand the derivation of the key Fokker- 
Planck equation given in Ref.0; we present some new re- 
sults for the evolving moments of the P2 (s, t) and correct 
some inconsequential errors in Ref . ! 



MODEL 

Our SOS model assigns an integer height h T to each 
point r on a square grid of dimensions L x x L y . Wc 
use periodic boundary conditions in the y direction. On 
our vicinal (001) simple cubic crystal, we create N close- 
packed [100] steps, with mean separation L = L x /N, via 
screw periodic boundary conditions in the x direction. In 
our simulations we take iV=5 in the initial simulations^ 
and iV=20 in later investigations. The energy of a config- 
uration is given by the standard absolute SOS prescrip- 
tion: 



H 



-E a ^ h r h r+ s 



(4) 



where S runs over the 4 nearest-neighbors of a site, and 
the factor 1/2 cancels the double-counting of bonds. 

In our SOS model, which has been described elsewhere 
fl2| . we use barriers determined by the standard long- 



standing simple rule[13|, [l4j, ll5[ of bond-counting: the 



barrier energy Eb is a diffusion barrier Ed plus a bond 
energy E a times the number of lateral nearest neighbors 
in the initial state. This number is 1 for an edge atom 
leaving a straight segment of step edge for the terrace, 
3 for a detaching atom that originally was part of this 
edge (leaving a notch or kink-antikink pair in the step), 
or 2 for a kink atom detaching, either to the step edge 
or the terrace. Processes that break 4 bonds, in partic- 
ular the removal of an atom from a flat terrace plane, 
are forbidden, as is any form of sublimation. No Ehrlich- 
Schwoebel barrier hinders atoms from crossing steps. We 
chose values 0.9 < E d < 1.1 eV and 0.3 < E a < 0.4 eV, 
using temperatures 520K <T< 580K. At these temper- 
atures we expect no significant finite-size effects in the y 
direction for the values of the mean terrace width L (in 
lattice spacings) that we use: 4 < L < 15. 

The width L y of the lattice should be greater than 
the "collision length" y co \\, the distance along y for 
a step to wander a distance L/2 in x. For a TSK 
model, estimates using a random- walk model give y co \\ = 
(L 2 /2)smh 2 (E k /2k B T) 461], where E k is the formation 
energy of a kink. At the temperatures and energies used 
in our simulation, y co n is of order 10 2 for L=6 and 10 3 
for L=15. E.g., for T=580K, E k = E a /2 = 0.175eV, 
and L=6, y C oii ~ 140. In almost all simulations reported 
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here, we use L y = 10 4 . While L y may often be larger 
than necessary, it allows for some self-averaging, decreas- 
ing the number of runs we need to carry out to get good 
statistics. 

In our rejection-free KMC, we separate all top-layer 
sites into 4 classes, those with z=0,l,2,3 nearest neighbors 
(NNs). (Those with z=4 are not allowed to move and are 
not considered when updating.) Typical realizations of 
these four classes are isolated adatoms, atoms protruding 
from a straight step edge, atoms at kink sites, and atoms 
at the edge of a step, respectively. We compute proba- 
bilities for each of the movable classes: Pi = fi/ J2i=o /*> 
where /j = Ni ■ cxp[— (Ed + i * E a )/kBT], and N is the 
number of sites with i NNs. (Of course, the 4 exponenti- 
ations are done once and for all at the beginning for each 
set of energies.) For each update we need four random 
numbers — ri, r%, r$, — uniformly distributed between 
and 1. We use r\ to pick which of the 4 movable classes 
will have the move. For the "winning" class, deter- 
mines which of the -/V; possible atoms will move. Then 
r3 determines in which of the 4 NN directions the atom 
moves. In this rejection-free scheme, we then decrease 
the height (the z value) of the initial position by one 
and increase the height of the chosen direction move 
from this initial site by one. This scheme can be (and 
has been, elsewhere) modified to allow for an Ehrlich- 
Schwoebel barrier. Finally is used to advance the 
clock in standard KMC fashion, similar to the n-fold 
way or bkl[I3i approach, using the prescription At = 
— hi(r4,)/R, where R is the total rate for a transition 
from the initial state [l(| • Explicitly, R = v$ Yn=o fi 
•-" , exp[-E d /k B T] J2i=o N 'i ■ exp[— i * E a /k B T], where we 
take the hopping frequency vq = 10 13 s _1 . 

We saved essentially every hundredth update; that in- 
terval corresponds to our unit of time, which is about 
1 sec. for the selected temperature and energies jl8| ■ This 
update interval is long enough so that the sum of the 
KMC update times varies insignificantly (±0.01%) but 
short enough to capture the behavior during the steep 
initial rise. 

In our model, the mass carriers are atoms rather than 
vacancies (or both). Since atoms with i=4 are frozen 
in our model, atom-vacancy pairs cannot form sponta- 
neously on a terrace. (More generally, this mechanism 
is highly improbable.) Mass carriers are thus created 
at step edges. If MC moves depend on the difference be- 
tween final and initial energies, as in Metropolis schemes, 
then there is equivalence between atom and vacancy cre- 
ation and transport. (If one goes beyond a strict SOS 
model and allows local relaxation, vacancies tend to be 
favored somewhat 0,[i^.) An atom quitting a step edge 
for the lower terrace costs 3E a if it leaves a straight step 
and 2E a if it leaves from a kink. At the upper side of a 
step, a vacancy can be spawned if a step-edge atom moves 
out one spacing onto the lower terrace (with the same en- 




FIG. 1: Three examples of fits using Eq. ((3]), used to extract 
r. Note that the data are very well fit in all three cases. The 
plotted a(t) is the standard deviation of the KMC data di- 
vided by the mean step spacing L (listed in lattice constants), 
and Ed and E a are the energy barriers for diffusion and for 
breaking a bond, respectively. Time t is essentially in seconds 
(see text). In all cases here and in later figures, q=2. Here 
AT=5. 



ergy cost as just given) and its inner neighbor happens to 
move in the same direction before the initial atom returns 
to its initial position. In kinetic Monte Carlo, however, 
rates are determined just by the difference between the 
barrier energy and the initial-state energy. This does not 
change the energy to produce an atom, but adds a cost 
of 3E a for the move of the second, inner-neighbor atom. 
Moreover, while the energy for an atom to hop along the 
terrace is Ed, for a vacancy it is Ed + 3E a . Indeed, we 
never observed the unlikely concerted process for vacancy 
creation in our simulations nor, for that matter, did we 
see any vacancies. The number of isolated atoms was 
also very small, with Nq being in single digits, and they 
moved very rapidly, rarely appearing in successive saved 
images. 

The freezing of i=4 processes marks a violation of de- 
tailed balance (since such a vacancy, if it existed, could be 
filled by a roving adatom); however, given the negligible 
occurrence of such vacancies in our simulations, the vio- 
lation should be insignificant. In some physical systems, 
motion of surface vacancies does evidently dominate mass 
2l|. Again, our goal in these calculations is 



transport 

not to account generally for experiments but to create 
a fully-controlled data set to see how well the dynamics 
can be described using our Fokker-Planck formalism. 



COMPUTED RESULTS 

We extract a characteristic time (or inverse rate) r 
from numerical data by fitting the dimensionless width 
using Eq. ([3|) , as illustrated in Fig. [TJ The fit is notably 
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(E d ,3E g )/k B T 

FIG. 2: Semilog plots of the relaxation time r (in sec.) vs. 
the diffusion barrier Ed (squares, upper line, red) or thrice the 
bond energy E a (triangles, blue), with the other held fixed, 
both in eV, with k B T = 0.05eV and iV=5. The numbers 
indicate the slopes; both are essentially unity. 

better than that found in the Metropolis/TSK study in 
Ref. 0- [However, the saturation value is notably higher 
than in Ref. 0, with the normalized standard deviation a 
(the value in the simulation divided by L) approaching 
^0.48, or a dimensionless variance of 0.24, rather than 
0.18 as found in Ref. H and anticipated from Eq. ((3]). 
This difference arises because the present algorithm al- 
lows steps to make contact along edge links rather than 
just at corners as in the usual fermion simulations. The 
variance of 0.18 is appropriate to "free fermions" with 
g=2. As we discuss in detail elsewhere the present 
algorithm leads to a smaller (and L-dcpendcnt) effective 
g as the steps come in contact more frequently, i.e. for 
smaller L and higher T (cf. Fig. [I]). For the present choice 
of parameters (i=6, k B T/ E d &l/20), the TWD has close 
to £>=1, for which the dimensionless variance is 0.27. This 
feature is inconsequential for the arguments in this pa- 
per.] 

We expect that the decay time exhibits Arrhenius be- 
havior: t cx cxp(Eb/ ksT). We investigate Et closely 
in the two traces of Fig. [2l We show typical runs at 
T = 580K, corresponding to k B T ps 1/20 eV. First, 
we ramped Ed, holding E a fixed at 0.35 eV (open 
squares, red). In the semi-log plot of reduced energies 
(energies/ ksT), we find a slope of 0.99 ± 0.02, indicat- 
ing that in the effective barrier, the multiplier of Ed/ksT, 
is essentially unity, as expected. In a second set of runs, 
we ramped E a , holding Ed fixed at 1.0 eV (open trian- 
gles, blue). Plotting now vs. 3E a /kBT, we find a slope 
0.94 ± 0.02, indicating that the effective energy barrier 
Et is Ed + 3E a . To corroborate this idea, we ramped T 
from 520 to 580K, fixing E a = 0.35eV and E d = l.OeV. 




FIG. 3: Demonstration of the robustness of the form of the 
evolving standard deviation, a(t) — ct(oo) [1 — exp(— t/r)] 1 ^ 2 , 
for five temperatures. In the inset, the five values of r 
by which the curves are rescaled are plotted vs. (Ed + 
3E a ) /ksT , showing their Arrhenius form. Here _Ed=l-0eV 
and _E a =0.35eV. In this and subsequent figures, iV=20. 

As illustrated in the inset of Fig. [3j we determine the 
fitted activation energy to be 2.03 ± 0.03, in excellent 
agreement with Ed + 3E a = 2.05 [eV]. Evidently the 
rate-determining process is the removal of a 3-bonded 
atom from a straight step, creating a pair of kinks (i.e., 
a kink and an antikink[23j) rather than the presumably 
more frequent process, with energy Ed + 2E„ , in which an 
atom leaves a kink position of a step [23j, [24| . (Of course, 
kink-antikink pairs also arise with a lower barrier when 
an atom from the terrace or from a kink site attaches to a 
step edge or splits off from a kink site. However, as mem- 
bers of class i = 1, such edge-atom structures are likely 
to be very short-lived.) The main part of Fig. [3] shows 
the standard deviation (oc TWD width) vs. time scaled 
by the relaxation time of each of the five temperatures. 

1/2 

Evidently the fit to a(t) = a(oo) [1 — exp(— t/r)] is 
robust. 

Initially the steps retreat as atoms are emitted. There 
is also an asymmetry in fluctuations from a straight step, 
since retreating moves involve higher barriers than ad- 
vancing fluctuations. Once the continuum picture be- 
comes applicable, the fluctuations appear to be symmet- 
ric, with typical configurations shown in Fig. 21 

To check consistency, we compare the intercepts of the 
linear fits in the two semilog plots, i.e., the prefactors 
of the exponential term in which the particular energy 
is ramped. In addition to the activation components 
there is the leading factor t = (w) 2 /4^o[2, EH) where 
we make the standard assignment for the hopping fre- 
quency, vq = 10 13 Hz. Since (w) = L = 6 in our simulations, 
tq, theoretically expected to be 9 x 10~ 13 s, is found in 
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FIG. 4: Typical step configurations during the evolution of 
initially straight steps in Fig. [3] The panels are 120 x 5000 
site portions extracted from the full 120 x 10000 net; there 
is considerable compression in the y direction. The panels 
range from early time to near saturation. Specifically, the 
ratios of the image time to r are: ~ 1/80, ~ 1/8, ~ 1/2, and 
somewhat over 3. From the step images alone, one would be 
hard pressed to distinguish the uphill direction, which is to 
the left. 



the simulations to be (8.71 ±0.25) x 10~ 13 s. In the ramp 
of Ed, the prcfactor is tq exp(3E a / ksT) , predicted to be 
1.19 x 10 _3 s. The value we find from the simulations 
is (1.2 ± 0.1) x 10~ 3 s, in excellent agreement. Similarly 
in the ramp of E a , the prefactor TQCxjp(Ed/ksT) is pre- 
dicted to be 4.366 x 10 _4 s and measured from the fit as 
(4.36 ±0.15) x 10~ 4 s. 

We also varied the system size L x , holding the num- 
ber of steps fixed, and thereby ramping (w). From the 
random- walk analogy, the prediction is that r oc (w) 2 . 
We find tolerable agreement, with a slope 18% below the 
expected value. We suspect that the reason behind the 
poorer agreement than for L=6 above originates in the 
L-dependence of the variance associated with the pecu- 
liar algorithm used in our simulations, which allows steps 
to touch [12] . 

Independently, another argument corroborates that 
the kink creation rate has an activation energy of E& + 
3E a : At equilibrium the creation and the annihilation 
rates are equal, so we compute the latter. Annihilation 
of kinks requires that an adatom diffuses to a step-edge 
notch — a kink-antikink pair, whose density is rik-ak ~ 
exp(-2E k /k B T) w exp (-2(E a /2)/k B T). Since the 
equilibrium adatom density is c cq = exp(— 2i? a /fcsT), 
the annihilation rate of kinks at a step edge is propor- 
tional to Dc cq n k _ ak ~ exp [-(Ed + 2E a + E a )/k B T] (cf. 
Ref. l26h . This in turn implies that the activation energy 




FIG. 5: Checks of dependencies on initial conditions and up- 
date moves. Evolution of the standard deviation a of the 
TWD for three initial configurations: straight steps (solid, 
black), "decimated" edge (dash-dotted, red), and crenelated 
(dotted, blue) edge. The 120 x 10,000 lattice has 20 steps, 
with L=6; we choose T=580K, E d =leV, and £' Q =0.4eV. For 
equilibrated non-interacting (free-fermion-like, g = 2) steps, 
a approaches 0.42 = aw, as expected. The smooth curve is 
a fit to an exponential approach to saturation. The smooth 
curve is <j(t) = cr(oo) [1 - exp(-i/74852)] 1/2 . Inset: Surface 
azimuthally misoriented by 0.0005 radians, forcing 5 kinks 
along the 10 4 -site steps. The three filled (red) circles repre- 
sent runs with an initial "perfect" configuration of 5 straight 
2000-site segments; the slope, indicated by the solid red line, 
is 1.0±0.1, in excellent agreement with the steeper line in 
Fig. 2. Thus, processes in which 3 bonds are broken gov- 
ern the scaling of the relaxation time r (given in sec. as in 
Fig. [1]). The open circles are from runs with 3-bonded atoms 
immobile; the corresponding slope is 0.67T0.03, so that now 
2-bonded atoms control the (much larger) t. 



for the kink creation rate is E c i + 3E a . 



DISCUSSION OF RESULTS 
Crucial role of kink creation 

The key energy in the relaxation time is that for de- 
taching 3-bonded atoms rather than kink atoms, a re- 
markable observation. Neither equilibrium nor growth 
processes involve 3-bonded atoms. At equilibrium, step 
fluctuations are controlled by the so-called step mobility, 
which isproportional to the emission rate of adatom from 
kinks 27J, which involves 2-bonded atoms. Hence, the re- 
laxation towards equilibrium should also be controlled by 
the step mobility. However, our results evidently contra- 
dict this notion. 

One possible explanation of our remarkable finding is 
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that kinks have to be formed first, requiring the extrac- 
tion of atoms from straight steps. (As noted earlier, the 
addition of an atom to a step edge also creates kink- 
antikink pairs, but the "tooth" is of class Nx, so very 
short-lived.). In that case, our initial configuration, in 
which steps are perfectly parallel and straight may be 
introducing a bias in the results. To investigate this pos- 
sibility, we considered two other initial states with equal 
numbers of kinks and antikinks, i.e. in which one pro- 
duces kinks by adding atoms to (or removing atoms from) 
a straight edge [28| . In the "decimated" case every tenth 
atom along a straight step is removed; in the other case 
every other atom is removed to create a "fully kinked" 
step, so that the edge resembles dentil molding or castle 
crenelations. In Fig. [5] we plot the resulting evolution of 
the standard deviation a of the TWD for the three cases. 
During the early-time rapid spreading of the initial sharp 
TWD, the slope increases with the number of initial 
kinks. In this regime, our model is not expected to apply 
(nor should any other ID, continuous model); indeed, 
Eq. ((3|) does not describe the steep initial part of the 
traces very well. After about 4x 10 4 MCS the curves are 
essentially indistinguishable within the noise level. The 
smooth curve is a one-parameter best fit of the data for 
the initially straight step by a(t)/a(oo) = 1 — cxp(— t/r). 
For this case we find r w 7.5T0 4 (in units that are essen- 
tially sec). 

We formulated the crenelated configuration because it 
creates at the outset a high density of atoms of class N\ : 
8 x 10~ 2 atoms/site, three orders of magnitude greater 
than the equilibrium density of 5 x 10~ 5 . These atoms 
quickly lead to a burst of adatoms that should quickly 
thermalize the step configurations. If it were only the 
supply of adatoms that limits equilibration, then for this 
scenario the subsequent creation of kinks would not be 
crucial. Evidently this is not so; even for the crenelated 
case, the relaxation kinetics are determined by the rate 
of creating kinks-antikinks pairs, with the usual energy 
barrier. 

To corroborate that kink creation is indeed the rate- 
limiting process, we computed the relaxation rate of a 
surface with steps azimuthally misoriented so as to create 
kinks via screw boundary conditions in the y direction. 
Specifically, in the initial state the in-plane misorienta- 
tion slope was set at 0.0005, so that geometry forces the 
existence of 5 kinks for L a =10,000. Keeping the diffusion 
barrier fixed at 1 eV, we varied E a . The results are shown 
in the inset of Fig. [5] as filled circles. We computed just 
three points, but clearly, essentially no difference is found 
with respect to the relaxation rate of straight [100] steps 
(the red line in Fig. [2]). The latter is drawn as a dashed 
line in the inset. The fitted slope to the data (times fceT) 
is 3.0 ± 0.3, fully consistent with 3-bonded ledge atoms 
being responsible for the rate-limiting process. We also 
checked that the relaxation rate is enormously slowed if 
3-bondcd atoms are kept immobile: Fitting the distri- 



bution width with Eq. ([3]), we ramped E a while hold- 
ing fixed Ed — 1 cV. The extracted relaxation times are 
shown in the inset of Fig. [5] as open circles. The reduced 
slope is 2.0 ±0.1, consistent with 2-bondcd kink atoms 
providing the rate- limiting process for the step motion in 
this case. The characteristic time is at least an order of 
magnitude larger than the previous case, which can be in- 
terpreted as due to the inability to create new kink sites, 
so that the number of sources for 2-bond escape of atoms 
to the straight segments of the step is limited to the ini- 
tial 5 kinks. Without the azimuthal misoricntation, this 
surface would be inert. Furthermore, the eventual width 
of the distribution, cr sat , is only about half the size of the 
3-bond case. Thus, at least over the course of our long 
runs, the surface is never able to equilibrate. 

We considered the number of N2 sites, typically kinks 
along steps. This quantity rose much more rapidly than 
the variance. Referring to the main plot of Fig. [5l N2 
achieves its saturation value by t k 3 x 10 4 . Thus, the 
process controlling r is not the initial formation of an 
adequate number of kinks but rather the maintenance 
of this number. For our chosen energies, about 1 in 30 
sites along a step was a kink, far higher than in our az- 
imuthally slightly-misoriented case. 

We also applied a similar analysis of the variance of the 
TWD to a vicinal (001) surface misoriented in along an 
azimuth rotated 45° so as to have zig-zag [110] steps. For 
such steps, every outer atom has i~2 lateral neighbors. 
Our analysis then shows that these 2-bond kink atoms 
produce the rate-limiting step, with a slope of 2 in the 
equivalent of the plot of r vs. E a in the inset of Fig. [5] 
This system has some idiosyncratic behavior due to the 
ease of creating fluctuations of steps from their mean 
configuration. Discussions of these subtleties would cloud 
the focus of this paper. Hence, we defer details to a future 
communication [22j. 

In short, we reach the striking conclusion that the equi- 
libration of a terrace width on a vicinal (001) simple cu- 
bic crystal with close-packed [100] steps (or steps not far 
from close-packed) and the fluctuations of the same ter- 
race width at equilibrium arc qualitatively different phe- 
nomena. The latter can take place with a constant num- 
ber of kinks, while the former requires creation of new 
kinks. Fluctuations from the equilibrium distribution, 
having the form of Eq. (2), will not lead to arbitrary ini- 
tial configurations such as a perfect cleaved crystal with 
straight, uniformly-spaced steps. 

Higher moments of the TWD 

So far, we have only considered the first two moments 
of the TWD. Several different distributions might ac- 
count for these two moments. As a further check that 
P2{s,t) describes the KMC data well, we study higher 
moments of the TWD in comparison with the analytical 
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0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 
f (10 5 ) 

FIG. 6: From top to bottom, the second (squares, red), third 
(circles, black), and fourth (triangles, blue) moments (with 
respect to the origin) of the evolving KMC-generated TWDs 
for initially straight configurations and the same parameters 
as in Fig. [5] for ease of comparison, the j moment fij (i) is 
divided by its equilibrium value /4f(oo). To make contact with 
the exponential approach of these moments to their saturation 
values, one must rescale the KMC time. For /i2 the rescaling 
factor, as in Fig. [5] is 7.5 TO 4 . For /13 and 114, the approach 
to saturation is progressively slower; the rescaling times are 
1.15-10 5 and 1.37-10 5 , 3/2 and 9/5 as large, respectively. 

expressions in Eqs. (|Al"6|) and (KT7\i . In Fig.[f5]we plot the 
second, third, and fourth moments (with respect to the 
origin) for initially-straight steps. (Since we wish to com- 
pare with P2(s, i), we divide the "raw" KMC j th moment 
by the j power of the mean spacing to determine /ij(i).) 
For each moment there is a steady rise (from unity) that 
approaches the equilibrium value fij(oo) exponentially. 
(Steps which arc initially decimated or crenelated behave 
similarly, though in somewhat "noisier" fashion.) For )i2, 
/13, and these saturation values agree well with the 
analytic results 3tt/8 ps 1.18, tt/2 ps 1.57, 15tt 2 /64 « 2.3, 
respectively, as which can be read off Eqs. (|A12|1 . (|A16|) . 
and (|A17[) . For ease of comparison, we plot /j, j (t) / fij (oo) 
in Fig. [HI so that each normalized moment approaches 
unity. The approach to saturation is evidently slower for 
successively higher moments. 

To make contact between the moments extracted from 
the KMC data and the analytic expressions in dimcn- 
sionlcss units arising from our Fokkcr-Planck analysis, we 
seek whether by rescaling KMC times by some tj leads to 
a good description of the data by the deduced moments. 
For [j,2 such a rescaling factor t2, with t2 = 7.5 • 10 4 , 
was already used in analyzing the data in Fig. [5] In 
other words, we adjust t2 so that /^(i/^V/^oo) from 
Eq. (|A12|) fits the data as well as possible, as illustrated 
in Fig. [6j For /j.^ and /X4, the approach to saturation is 
progressively slower; the rescaling times £3 and £4, simi- 
larly obtained, are 1.15-10 5 and 1.37-10 5 , 3/2 and 9/5 as 




FIG. 7: Analytic results for the moments of the evolving 
TWD predicted by the (ID continuum) Fokker-Planck the- 
ory: Plot of the effective exponential decay time of the dif- 
ference between the evolving second (solid, red), third (long- 
short-dashed, black), and fourth (dashed, blue) moments of 
the solution to Eq. JT}, P e= 2(s,f), given explicitly in the ap- 
pendix in Eqs. (|A12f) . (|A16|I . and (|A17|I . and their steady- 
state, asymptotic values associated with Eq. ((3]). As in the 
KMC data in Fig. [6] the decay is significantly slower for the 
higher moments. The thicker set of curves corrects for the 
modest deficiency of the the first moment (J.i(t); Hi(t) is de- 
picted in the inset (upper left, italicized axes labels) and given 
analytically in Eq. l|A13[l . See text for details. 



large, respectively. 

The analytic expressions for the four moments, given in 
Eqs. (|A12j) . (|A13|) . (|A16|) . and (|A17j) . all approach satu- 
ration asymptotically from below like cxp(— t). However, 
in the temporal regime corresponding to the KMC simu- 
lations, higher-order terms cause the evident exponential- 
like approach to depend on t/r rather than simply t, 
where r is some effective time constant of order unity. 

To determine Tj we consider in Fig. [7] the evolution of 
— f/ln[l — /j,j(t)/ Hj(oo)] for each moment (j=2,3,4). To 
the degree that this trace is horizontal, the Ansatz is 
appropriate. The thin curves in Fig. [7]show that this as- 
sumption becomes progressively better as time advances. 
While all three curves eventually converge to unity, in the 
time regime under consideration the three moments have 
significantly different time constants, with the higher mo- 
ments having progressively larger magnitudes, consistent 
with the KMC findings displayed in Fig. [6l 

The inset of Fig. [7] displays the first moment of fJ,i(t). 
As described in the appendix, it is a couple percent 
smaller than the proper value of unity in the region 
around t « 1. (This is clearly a deficiency only of the 
analytic results. The KMC data have Hi(t) = 1 by con- 
struction.) Analogous to the transformation of the TWD 
from a function of w to P(s, t), where s = w/(w), we can 
rewrite the distribution in terms of s/ fix (t) and show that 
then the "corrected" moments //™ rr (i) = fJ.j{t)/ fJ-\(t)- 
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Obviously, /if T (t) = 1. The higher moments yu§ orr (t), 
/i§ orr (t), and /Z4° rr (f) arc displayed as the thick curves 
in Fig. [7l These curves flatten considerably sooner than 
the thin curves, and to a value ~ 1/2, reminiscent of 
Eq. (|A15[) . For these curves, (r 3 — t 2 )/t 2 is somewhat 
over 0.1 while (Y4 — t 2 )/t 2 is about twice as large, cap- 
turing the trend of the numerical data. If we use the ana- 
lytic expressions for the corrected moments as the rescal- 
ing factors, then the time rescaling factors are 1.75 -10 5 , 
1.9-10 5 , and 2.0-10 5 , respectively. Then (t 3 -t 2 )/t 2 « 0.09 
and (Y4 — t 2 )/t 2 ss 0.14, closer to the KMC values. 

Accounting for the skewness and the kurtosis poses a 
more difficult test for our kinetic Monte Carlo simulation 
of our model. Due to the interplay of several moments 
in computing these statistics, our numerical data is not 
adequate to test these predicted behaviors meaningfully. 
Such analyses are arguably the most stringent tests, in 
which we seek differences from random-walk, Gaussian 
behavior; they corroborate the conclusion above that the 
surface never fully equilibrates. Far more extensive com- 
putations might clarify this matter, but are beyond the 
scope of our present analysis. 

COMPARISON WITH GRUBER-MULLINS 

To help place our approach in context, we compare it 
with previous approximations in the literature based on 
the fermion description of the fluctuating steps. The cele- 
brated Gruber-Mullins (GM) approximation consid- 
ers a fluctuating step between two fixed neighbors treated 
as rigid boundaries. In fermion language, the step is the 
ID trajectory of a quantum particle confined to the seg- 
ment (0,2(uj)) by an infinite potential. Two cases are 
easily treated: For non-interacting steps, the fluctuating 
step is then the trajectory of a free fermion, and is equiv- 
alent to a classical particle performing a random walk in 
a potential of the form[3Cl] 

V{x) = -2ffi[sin(7ra/2(w))]. (5) 

Then s = w/(w) obeys the Langevin equation 

S (w) 2 tan(7rs/2) V ' [ > 

This approximation preserves the logarithmic behavior of 
the repulsive potential at short range; even the amplitude 
is correct: (2b e s — g/ s)\ e=2 — (8s/tt — 2/s) in Eq. (1) is 
replaced by — it cot(7rs/2), nearly the same for ,s < 0.7. 
However, the GM potential has bogus symmetry about 
(w), truncating the long-range tail of P(s). 

For strongly interacting steps, the fluctuating step feels 
an (approximately) quadratic confining potential, and 
the TWD distribution is predicted to be Gaussian. [29| . 
In our formalism the replacement is now (s— 1)/<Tq, where 
<Jq is the variance of the Gaussian TWD (so half the vari- 
ance of the associated ground-state wavefunction, which 



we used to construct the FP potential [3(j.) This replace- 
ment should be compared with (2b e s—g/s) ~ (g+^)s— g/s 
[31I ] . In the GM approximation, a^ 2 = [I2g(g — 2)] 1 ^ 2 , 
with (12) 1 / 2 « 3.5 replaced by (2tt 4 /15) 1 / 2 w 3.6 if the in- 
teractions with all steps rather than just the two bound- 
ing steps are considered [3lj |. An improved approxima- 
tion ("modified Grenoble") gives (Tq 2 r* 2.1g [31|. Our 
expression for the FP potential now differs at small s 
from that derived for these approximations because the 
Gaussians actually extend (unphysically, albeit with in- 
significant amplitude) to negative values of s. Our ap- 
proach is globally superior to the celebrated GM approx- 
imation (as well as to the usual alternatives [lj), both 
quantitatively and qualitatively, for all physical values 
of the step-step interaction strength [3^. Moreover, our 
FP equation ((T|) is fully soluble, so that the TWD can be 
obtained analytically as a function of time. 

SUMMARY 

In summary, we have shown that the relaxation time 
of the variance of the solution of our Fokker-Planck equa- 
tion for step relaxation on a vicinal surface can be fit to 
the comparable variance in a kinetic Monte Carlo simula- 
tion of the standard simple model of atomistic processes 
at surfaces. This time has Arrhcnius behavior that is 
related to microscopic processes, substantiating that this 
FP approach can offer useful physical insight into the evo- 
lution of complex surface structures toward equilibrium. 
Thus, once the continuum formalism becomes appropri- 
ate, the FP time in some sense tracks actual time in our 
model of an evolving physical system of steps with no 
energetic repulsion. The formalism also readily allows 
such repulsions, inviting future simulations to test how 
well the Fokker-Planck formalism describes such systems. 
Since the steps communicate from the outset, the contin- 
uum formalism might apply sooner. For the situation we 
have considered, we have argued in several ways that the 
rate-determining process in step relaxation is the creation 
of kink-antikink pairs. We have also examined higher 
moments of the distribution, both analytically and with 
simulations. While we make no pretense that our ap- 
proach is either exact or a formal theory, we have shown 
that it can be a fruitful way to treat relaxation of steps 
on surfaces. Many avenues of extension are possible. 
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Appendix: Derivation of Eq. ([T]) and Some 
Consequent Results 

In this appendix we expand the derivation of the 
Fokker-Planck equation given in Ref. 1 as well as cor- 
recting some algebra oversights in intermediate steps pre- 
sented there. We also present some simpler expressions 
for quantities of interest that arise for the case of non- 
interacting [energetically] steps {g = 2) investigated in 
the reported computations. 

As in Ref. 0, we begin with the correspondence found 
by Dyson between RMT and his Coulomb gas model § : 
iV classical particles on a line, interacting with a loga- 
rithmic potential, and confined by an overall harmonic 
potential. Dyson's model helps our understanding of the 
fluctuation properties of the spectrum of complex con- 
served systems. This model can be generalized to the 
dynamic Brownian motion model, in which the N parti- 
cles are subject, besides the mutual Coulomb repulsions, 
to dissipative forces [33j ]. The particle positions Xi then 
obey Langevin equations, 



E 

1^3 



(Al) 



where 77 is a delta-correlated white noise and g (oc g) is 
the "charge" of each particle. The probability of finding 
the particles at the positions {x n } at time t is the solution 
of the multidimensional FPE 



dP({x n },t) 
dt 



T- 



_d_ 

dx, 



P{{x n },t)+~/XiP({x n },t) 







rrf dxi 



-P({Xn},t) 



(A2) 



In the ID case, 7 -1 would essentially be the variance 
of the stationary distribution. Narayan and Shastry 
Q showed that the CS model is equivalent to Dyson's 
Brownian motion model, in the sense that the solu- 
tion of the FPE (|A2|) may be written as P({x n },t) = 
ip({x n },t)ipo({x n },t), where ip({x n },t) is the solution 
of a Schrodingcr equation with imaginary time, derived 
from the CS Hamiltonian. The deterministic force of 
Eq. JAQ) 



F(x m ) = -jx m - ^2 
so that 



E 



k>m q<m H 



(A3) 



E 



E 



X m +l 



_^\ X k Xm+l){Xk X m ) q<m ( x m+l X q )(x m X q ) 

Our goal is to find the distribution of widths w. Mind- 
ful of the Gruber-Mullins approach [2{|, we construct 
a single- "particle," mean-field approximation in which 
the dynamical variable is the nearest-neighbor distance 
w m = x rn+ i — x m . To decouple the force on w m from the 
other particles, we assume — in the spirit of GM — that 
the denominators (xh — x m +i){xk — x m ) in Eq. (|A4[) are 
replaced by their mean values, the average being taken 
in the stationary state: 

((x k -x m+1 )(x k -x m )) st = (w 2 ) st (k-m-l)(k-m), (A5) 



Each of the two sums in Eq. 
taking the form 



(x 



m+l 







' N 1 

E^TT)^ 



4| then simplifies greatly, 



N 



N + 1 N- 



(A6) 

Hence, the interaction of a particle pair with all other 
particles acts on average as a harmonic potential, increas- 
ing the "spring constant" of the external confining poten- 
tial. We arrive at a single-particle Langevin equation for 
the terrace width w: 



dii' 
~dt 



= -2 



7 Q 
2 (wi) st 



2Tr]. (A7) 



Our goal is to convert Eq. (|A7j) into a FPE for which 
Eq. ((3]) is a steady-state solution. We change to dimcn- 
sionless variables s = w/(w} s t an( l t = Tt/(w) 2 st . Treat- 
ing 7 as a self-consistency parameter and recognizing 
g = gT/2, we set 7 = T/(w 2 ) st - Then the coefficient 
in parentheses in Eq. (|A7|) becomes 



(A8) 



using the second moment of P e (s) [(s 2 ) = (g+ l)/(26 e )]. 
Furthermore, if (r](t)r](t')) = 5(t - t'), then rj(t) = 
y/2/T(w) 3t r)(t) satisfies (fj(t)fj(t')) = 5(t-t'). With these 
results, we recast Eq. (| A7|) into the Langevin equation 



di 



2b e s 



(A9) 



and thence the sought-after FPE given in Eq. {]]). 

To solve Eq. |T]) we must specify the initial distribution 
in so. For an initial (at t = 0) sharp distribution 8{s— 1), 
corresponding to a perfectly cleaved crystal, the solution 
is essentially written down by Montroll and West:[3^.[35j 



F(x m +i)-F(x m ) = -^f(x m+1 -x m ) - g 



(A4) 



P(s,t) = 2~b e s- 



(e-W 

e 4 Ig-1 



(2b e se-i^ 



(A10) 
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where b g = b g /(l — e - '). In the limit of long times, we 
showed in Ref. that, as t increases, this P(s, t) ap- 
proaches P e {s) of Eq. ©. 

For the particular case Q=2, b g becomes 4/[7r(l— e _t )], 
while Ii(z) = y/2/(irz) sinh(z). Then Eq. (fXlO]l simpli- 
fies to 



22 e * s 



7rsinh 2 (|) 



sinh 



' (4/tt)s 
v sinh(|) 





' 4(s 2 +c- f ~\ 


jexp 


7r(l-e-*f 



y.i) 



In experiments, P(s) is generally characterized just by 
its variance c 2 = fi2 — A* 2 , which can be calculated from 
its second and first moments, /Z2 and a*i, respectively: 



,-. 37T 

AW) = -^ u t 



-t / 3?r 
e = 1 



1 - e" 



1 (A12) 



In any case, the Arrhenius behavior of the characteris- 
tic time of the exponential will not be affected by such 
modest changes in the prefactor. 

If the Fokker-Planck description step relaxation is ro- 
bust, then the higher moments of P(s, t) should also char- 
acterize those moments extracted from the KMC data, as 
displayed in the text in Fig. [5J Thus, we present explicit 
analytic formulas for the third and fourth moments: 



M3(*~) 



1 _t 1/2 

2° U i 



5tt 3/2 
— u~ cxp 
16 * 



4/tt 



-e" 2t + 3ii,~e"' 



1 1 



-t/T 3 



(Al6a) 



Mi(f) = £ 



u~ exp 



-4/tt 



1- 



7T 



T(t) 



where, 
which 
tially 



(A13) 

for brevity, we take uj = 1 — exp(— i), 
obviously approaches unity exponcn- 
from below. Furthermore, we write 

T(t) = (7r/4)exp(i/2)erf{2/[7r(exp(t) - 1)] 1/2 ), where 
erf is the error function: [361] T(i) also approaches unity 
exponentially, but from above, after rising initially from 
7r/4 to about 1.01. 

Scrutiny of Eq. (|A13[) reveals that each of the two sum- 
mands in the square brackets approaches 1 for large t. 
As t approaches 0, the first summand vanishes while the 
second rises to 2. Thus, Hi(t) has the expected value for 
vanishing and large i, as illustrated in the inset of Fig. [7] 
There is, however, an initial rapid drop, reaching a min- 
imum of about 0.9745 around £=0.582, and then rising 
smoothly, reaching 0.99 by £=1.95, 0.995 by A=2.685, and 
0.999 by i=4.33. The small deviation from unity is pre- 
sumably due to the approximations in using Eq. (|A5j) to 
reach Eq. (|A7[) , which apparently break the symmetry of 
the fluctuations of the steps (m and m+1) bounding w m 

To the extent that this deviation is negligible [and in 
any case for qualitative purposes], we get 



a 2 (f)| 



mi=i : 



(A14) 



If we numerically evaluate cr 2 (t) using Eqs. (|A12|) and 
(IA13|) , we find a similar expression but with a more rapid 
rise to the equilibrium result; remarkably, it is well ap- 
proximated by 



a 2 (t)=a 2 ,(l- 



(A15) 



reminiscent of the solution of the Fokker-Planck equa- 
tion for a Brownian particle in a quadratic potential [381 ] . 



/i 4 (1) = 



157T 2 2 57T r 

-u- t + — e u~ t + e 



04 



15 
— 1 
64 



1 1 



-i/T t 



(A17) 
(Al7a) 



The approximate expressions for ^{i) and m\t) in 
Eqns. (|Al6aj) and (jAl7a| are written in the form of 
the exact result for ^(i) in Eq. (|A12jl . By setting 
T.3 = T4 = 1 one obtains a mediocre approximation which 
underestimates fis(t) and Hiit) by as much as 6% and 
15%, respectively. A far better accounting is obtained by 
taking T 3 « 0.79 and T4 m 0.76; the best values of these 
time constants depends weakly on the temporal range 
over which one seeks to optimize the agreement. The 
approximate expressions then underestimate the actual 
Hj(t) (by at most 2% and 3%) up to t w 3/2 and then 
overestimate it (by at most i% and 1%), respectively. 
Thus, H3(t) and Hi{t) can be well described by curves 
starting rising smoothly from unity and decaying expo- 
nentially toward their long-time limit, but with values of 
Tj that are smaller than unity. Finally, note that the ap- 
proximate expressions based on Eqns. (|A16aj) and (|A17ap 
are not used in the analysis of the moments in Subsection 
; hence, the values of the Tj play no role there. 

From these results the skewness can be expressed ana- 
lytically but has an unwieldy form. However, it is semi- 
quantitatively described by 0.4857 tanh(F) (i.e. to within 
±4% for t > 0.46 and within a percent for t > 1.9). In 
other words, the skewness rises smoothly and monotoni- 
cally from initially to the equilibrium value. Since for 
large t, tanh(A) ~ 1 — 2cxp(— 2t) we find the same ap- 
proach to saturation as for the variance in Eq. (|A15|) . The 
kurtosis begins at 3 but dips (to about 2.87 near t=0.5) 
before rising to its equilibrium value of 3.1082. The ap- 
proach to this asymptotic value is well approximated by 
3.11 (l-0.58e" 2t ). 
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